/// \defgroup newmat Newmat matrix manipulation library
///@{

/// \file newmat.h
/// Definition file for matrix library.

// Copyright (C) 2004: R B Davies

#ifndef NEWMAT_LIB
#define NEWMAT_LIB 0

#include "include.h"

#include "myexcept.h"

#include "gexport/ExportLibGnut.h"

#include "stdint.h"


#ifdef use_namespace
namespace NEWMAT { using namespace RBD_COMMON; }
namespace RBD_LIBRARIES { using namespace NEWMAT; }
namespace NEWMAT {
#endif

	//#define DO_REPORT                     // to activate REPORT

#ifdef NO_LONG_NAMES
#define UpperTriangularMatrix UTMatrix
#define LowerTriangularMatrix LTMatrix
#define SymmetricMatrix SMatrix
#define DiagonalMatrix DMatrix
#define BandMatrix BMatrix
#define UpperBandMatrix UBMatrix
#define LowerBandMatrix LBMatrix
#define SymmetricBandMatrix SBMatrix
#define BandLUMatrix BLUMatrix
#endif

// ************************** general utilities ****************************/

	class GeneralMatrix;                            // defined later
	class BaseMatrix;                               // defined later
	class MatrixInput;                              // defined later

	void MatrixErrorNoSpace(const void*);           ///< test for allocation fails

	/// Return from LogDeterminant function.
	/// Members are the log of the absolute value and the sign (+1, -1 or 0)
	class LibGnut_LIBRARY_EXPORT LogAndSign
	{
		Real log_val;
		int sign_val;
	public:
		LogAndSign() { log_val = 0.0; sign_val = 1; }
		LogAndSign(Real);
		void operator*=(Real);                       ///< multiply by a real
		void pow_eq(int k);                          ///< raise to power of k
		void PowEq(int k) { pow_eq(k); }
		void ChangeSign() { sign_val = -sign_val; }
		void change_sign() { sign_val = -sign_val; } ///< change sign
		Real LogValue() const { return log_val; }
		Real log_value() const { return log_val; }   ///< log of the absolute value
		int Sign() const { return sign_val; }
		int sign() const { return sign_val; }        ///< sign of the value
		Real value() const;                          ///< the value
		Real Value() const { return value(); }
		FREE_CHECK(LogAndSign)
	};

	// the following class is for counting the number of times a piece of code
	// is executed. It is used for locating any code not executed by test
	// routines. Use turbo GREP locate all places this code is called and
	// check which ones are not accessed.
	// Somewhat implementation dependent as it relies on "cout" still being
	// present when ExeCounter objects are destructed.

#ifdef DO_REPORT

	class ExeCounter
	{
		int line;                                    // code line number
		int fileid;                                  // file identifier
		long nexe;                                   // number of executions
		static int nreports;                         // number of reports
	public:
		ExeCounter(int, int);
		void operator++() { nexe++; }
		~ExeCounter();                               // prints out reports
	};

#endif


	// ************************** class MatrixType *****************************

	/// Find the type of a matrix resulting from matrix operations.
	/// Also identify what conversions are permissible.
	/// This class must be updated when new matrix types are added.

	class LibGnut_LIBRARY_EXPORT MatrixType
	{
	public:
		enum Attribute {
			Valid = 1,
			Diagonal = 2,             // order of these is important
			Symmetric = 4,
			Band = 8,
			Lower = 16,
			Upper = 32,
			Square = 64,
			Skew = 128,
			LUDeco = 256,
			Ones = 512
		};

		enum {
			US = 0,
			UT = Valid + Upper + Square,
			LT = Valid + Lower + Square,
			Rt = Valid,
			Sq = Valid + Square,
			Sm = Valid + Symmetric + Square,
			Sk = Valid + Skew + Square,
			Dg = Valid + Diagonal + Band + Lower + Upper + Symmetric
			+ Square,
			Id = Valid + Diagonal + Band + Lower + Upper + Symmetric
			+ Square + Ones,
			RV = Valid,     //   do not separate out
			CV = Valid,     //   vectors
			BM = Valid + Band + Square,
			UB = Valid + Band + Upper + Square,
			LB = Valid + Band + Lower + Square,
			SB = Valid + Band + Symmetric + Square,
			KB = Valid + Band + Skew + Square,
			Ct = Valid + LUDeco + Square,
			BC = Valid + Band + LUDeco + Square,
			Mask = ~Square
		};


		static int nTypes() { return 13; }          // number of different types
								// exclude Ct, US, BC
	public:
		int attribute;
		bool DataLossOK;                            // true if data loss is OK when
													// this represents a destination
	public:
		MatrixType() : DataLossOK(false) {}
		MatrixType(int i) : attribute(i), DataLossOK(false) {}
		MatrixType(int i, bool dlok) : attribute(i), DataLossOK(dlok) {}
		MatrixType(const MatrixType& mt)
			: attribute(mt.attribute), DataLossOK(mt.DataLossOK) {}
		void operator=(const MatrixType& mt)
		{
			attribute = mt.attribute; DataLossOK = mt.DataLossOK;
		}
		void SetDataLossOK() { DataLossOK = true; }
		int operator+() const { return attribute; }
		MatrixType operator+(MatrixType mt) const
		{
			return MatrixType(attribute & mt.attribute);
		}
		MatrixType operator*(const MatrixType&) const;
		MatrixType SP(const MatrixType&) const;
		MatrixType KP(const MatrixType&) const;
		MatrixType operator|(const MatrixType& mt) const
		{
			return MatrixType(attribute & mt.attribute & Valid);
		}
		MatrixType operator&(const MatrixType& mt) const
		{
			return MatrixType(attribute & mt.attribute & Valid);
		}
		bool operator>=(MatrixType mt) const
		{
			return (attribute & ~mt.attribute & Mask) == 0;
		}
		bool operator<(MatrixType mt) const         // for MS Visual C++ 4
		{
			return (attribute & ~mt.attribute & Mask) != 0;
		}
		bool operator==(MatrixType t) const
		{
			return (attribute == t.attribute);
		}
		bool operator!=(MatrixType t) const
		{
			return (attribute != t.attribute);
		}
		bool operator!() const { return (attribute & Valid) == 0; }
		MatrixType i() const;                       ///< type of inverse
		MatrixType t() const;                       ///< type of transpose
		MatrixType AddEqualEl() const               ///< add constant to matrix
		{
			return MatrixType(attribute & (Valid + Symmetric + Square));
		}
		MatrixType MultRHS() const;                 ///< type for rhs of multiply
		MatrixType sub() const                      ///< type of submatrix
		{
			return MatrixType(attribute & Valid);
		}
		MatrixType ssub() const                     ///< type of sym submatrix
		{
			return MatrixType(attribute);
		}        // not for selection matrix
		GeneralMatrix* New() const;                 ///< new matrix of given type
		GeneralMatrix* New(int, int, BaseMatrix*) const;
		///< new matrix of given type
		const char* value() const;                  ///< type as char string
		const char* Value() const { return value(); }
		LibGnut_LIBRARY_EXPORT friend bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
		LibGnut_LIBRARY_EXPORT friend bool Compare(const MatrixType&, MatrixType&);
		///< compare and check conversion
		bool is_band() const { return (attribute & Band) != 0; }
		bool is_diagonal() const { return (attribute & Diagonal) != 0; }
		bool is_symmetric() const { return (attribute & Symmetric) != 0; }
		bool CannotConvert() const { return (attribute & LUDeco) != 0; }
		// used by operator== 
		FREE_CHECK(MatrixType)
	};


	// *********************** class MatrixBandWidth ***********************/

	///Upper and lower bandwidths of a matrix.
	///That is number of diagonals strictly above or below main diagonal,
	///e.g. diagonal matrix has 0 upper and lower bandwiths.
	///-1 means the matrix may have the maximum bandwidth. 
	class LibGnut_LIBRARY_EXPORT MatrixBandWidth
	{
	public:
		int lower_val;
		int upper_val;
		MatrixBandWidth(const int l, const int u) : lower_val(l), upper_val(u) {}
		MatrixBandWidth(const int i) : lower_val(i), upper_val(i) {}
		MatrixBandWidth operator+(const MatrixBandWidth&) const;
		MatrixBandWidth operator*(const MatrixBandWidth&) const;
		MatrixBandWidth minimum(const MatrixBandWidth&) const;
		MatrixBandWidth t() const { return MatrixBandWidth(upper_val, lower_val); }
		bool operator==(const MatrixBandWidth& bw) const
		{
			return (lower_val == bw.lower_val) && (upper_val == bw.upper_val);
		}
		bool operator!=(const MatrixBandWidth& bw) const { return !operator==(bw); }
		int Upper() const { return upper_val; }
		int upper() const { return upper_val; }
		int Lower() const { return lower_val; }
		int lower() const { return lower_val; }
		FREE_CHECK(MatrixBandWidth)
	};


	// ********************* Array length specifier ************************/

	/// This class is used to avoid constructors such as
	/// ColumnVector(int) being used for conversions.
	/// Eventually this should be replaced by the use of the keyword "explicit".

	class LibGnut_LIBRARY_EXPORT ArrayLengthSpecifier
	{
		int v;
	public:
		int Value() const { return v; }
		int value() const { return v; }
		ArrayLengthSpecifier(int l) : v(l) {}
	};

	// ************************* Matrix routines ***************************/


	class MatrixRowCol;                             // defined later
	class MatrixRow;
	class MatrixCol;
	class MatrixColX;

	class GeneralMatrix;                            // defined later
	class AddedMatrix;
	class MultipliedMatrix;
	class SubtractedMatrix;
	class SPMatrix;
	class KPMatrix;
	class ConcatenatedMatrix;
	class StackedMatrix;
	class SolvedMatrix;
	class ShiftedMatrix;
	class NegShiftedMatrix;
	class ScaledMatrix;
	class TransposedMatrix;
	class ReversedMatrix;
	class NegatedMatrix;
	class InvertedMatrix;
	class RowedMatrix;
	class ColedMatrix;
	class DiagedMatrix;
	class MatedMatrix;
	class GetSubMatrix;
	class ReturnMatrix;
	class Matrix;
	class SquareMatrix;
	class nricMatrix;
	class RowVector;
	class ColumnVector;
	class SymmetricMatrix;
	class UpperTriangularMatrix;
	class LowerTriangularMatrix;
	class DiagonalMatrix;
	class CroutMatrix;
	class BandMatrix;
	class LowerBandMatrix;
	class UpperBandMatrix;
	class SymmetricBandMatrix;
	class LinearEquationSolver;
	class GenericMatrix;


#define MatrixTypeUnSp 0
	//static MatrixType MatrixTypeUnSp(MatrixType::US);
	//						// AT&T needs this

	/// Base of the matrix classes.
	class LibGnut_LIBRARY_EXPORT BaseMatrix : public Janitor
	{
	protected:
		virtual int search(const BaseMatrix*) const = 0;
		// count number of times matrix is referred to
	public:
		virtual GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp) = 0;
		// evaluate temporary
// for old version of G++
//   virtual GeneralMatrix* Evaluate(MatrixType mt) = 0;
//   GeneralMatrix* Evaluate() { return Evaluate(MatrixTypeUnSp); }
		AddedMatrix operator+(const BaseMatrix&) const;    // results of operations
		MultipliedMatrix operator*(const BaseMatrix&) const;
		SubtractedMatrix operator-(const BaseMatrix&) const;
		ConcatenatedMatrix operator|(const BaseMatrix&) const;
		StackedMatrix operator&(const BaseMatrix&) const;
		ShiftedMatrix operator+(Real) const;
		ScaledMatrix operator*(Real) const;
		ScaledMatrix operator/(Real) const;
		ShiftedMatrix operator-(Real) const;
		TransposedMatrix t() const;
		//   TransposedMatrix t;
		NegatedMatrix operator-() const;                   // change sign of elements
		ReversedMatrix reverse() const;
		ReversedMatrix Reverse() const;
		InvertedMatrix i() const;
		//   InvertedMatrix i;
		RowedMatrix as_row() const;
		RowedMatrix AsRow() const;
		ColedMatrix as_column() const;
		ColedMatrix AsColumn() const;
		DiagedMatrix as_diagonal() const;
		DiagedMatrix AsDiagonal() const;
		MatedMatrix as_matrix(int, int) const;
		MatedMatrix AsMatrix(int m, int n) const;
		GetSubMatrix submatrix(int, int, int, int) const;
		GetSubMatrix SubMatrix(int fr, int lr, int fc, int lc) const;
		GetSubMatrix sym_submatrix(int, int) const;
		GetSubMatrix SymSubMatrix(int f, int l) const;
		GetSubMatrix row(int) const;
		GetSubMatrix rows(int, int) const;
		GetSubMatrix column(int) const;
		GetSubMatrix columns(int, int) const;
		GetSubMatrix Row(int f) const;
		GetSubMatrix Rows(int f, int l) const;
		GetSubMatrix Column(int f) const;
		GetSubMatrix Columns(int f, int l) const;
		Real as_scalar() const;                      // conversion of 1 x 1 matrix
		Real AsScalar() const;
		virtual LogAndSign log_determinant() const;
		LogAndSign LogDeterminant() const { return log_determinant(); }
		Real determinant() const;
		Real Determinant() const { return determinant(); }
		virtual Real sum_square() const;
		Real SumSquare() const { return sum_square(); }
		Real norm_Frobenius() const;
		Real norm_frobenius() const { return norm_Frobenius(); }
		Real NormFrobenius() const { return norm_Frobenius(); }
		virtual Real sum_absolute_value() const;
		Real SumAbsoluteValue() const { return sum_absolute_value(); }
		virtual Real sum() const;
		virtual Real Sum() const { return sum(); }
		virtual Real maximum_absolute_value() const;
		Real MaximumAbsoluteValue() const { return maximum_absolute_value(); }
		virtual Real maximum_absolute_value1(int& i) const;
		Real MaximumAbsoluteValue1(int& i) const
		{
			return maximum_absolute_value1(i);
		}
		virtual Real maximum_absolute_value2(int& i, int& j) const;
		Real MaximumAbsoluteValue2(int& i, int& j) const
		{
			return maximum_absolute_value2(i, j);
		}
		virtual Real minimum_absolute_value() const;
		Real MinimumAbsoluteValue() const { return minimum_absolute_value(); }
		virtual Real minimum_absolute_value1(int& i) const;
		Real MinimumAbsoluteValue1(int& i) const
		{
			return minimum_absolute_value1(i);
		}
		virtual Real minimum_absolute_value2(int& i, int& j) const;
		Real MinimumAbsoluteValue2(int& i, int& j) const
		{
			return minimum_absolute_value2(i, j);
		}
		virtual Real maximum() const;
		Real Maximum() const { return maximum(); }
		virtual Real maximum1(int& i) const;
		Real Maximum1(int& i) const { return maximum1(i); }
		virtual Real maximum2(int& i, int& j) const;
		Real Maximum2(int& i, int& j) const { return maximum2(i, j); }
		virtual Real minimum() const;
		Real Minimum() const { return minimum(); }
		virtual Real minimum1(int& i) const;
		Real Minimum1(int& i) const { return minimum1(i); }
		virtual Real minimum2(int& i, int& j) const;
		Real Minimum2(int& i, int& j) const { return minimum2(i, j); }
		virtual Real trace() const;
		Real Trace() const { return trace(); }
		Real norm1() const;
		Real Norm1() const { return norm1(); }
		Real norm_infinity() const;
		Real NormInfinity() const { return norm_infinity(); }
		virtual MatrixBandWidth bandwidth() const;  // bandwidths of band matrix
		virtual MatrixBandWidth BandWidth() const { return bandwidth(); }
		void IEQND() const;                         // called by ineq. ops
		ReturnMatrix sum_square_columns() const;
		ReturnMatrix sum_square_rows() const;
		ReturnMatrix sum_columns() const;
		ReturnMatrix sum_rows() const;
		virtual void cleanup() {}
		void CleanUp() { cleanup(); }

		//   virtual ReturnMatrix Reverse() const;       // reverse order of elements
		//protected:
		//   BaseMatrix() : t(this), i(this) {}

		friend class GeneralMatrix;
		friend class Matrix;
		friend class SquareMatrix;
		friend class nricMatrix;
		friend class RowVector;
		friend class ColumnVector;
		friend class SymmetricMatrix;
		friend class UpperTriangularMatrix;
		friend class LowerTriangularMatrix;
		friend class DiagonalMatrix;
		friend class CroutMatrix;
		friend class BandMatrix;
		friend class LowerBandMatrix;
		friend class UpperBandMatrix;
		friend class SymmetricBandMatrix;
		friend class AddedMatrix;
		friend class MultipliedMatrix;
		friend class SubtractedMatrix;
		friend class SPMatrix;
		friend class KPMatrix;
		friend class ConcatenatedMatrix;
		friend class StackedMatrix;
		friend class SolvedMatrix;
		friend class ShiftedMatrix;
		friend class NegShiftedMatrix;
		friend class ScaledMatrix;
		friend class TransposedMatrix;
		friend class ReversedMatrix;
		friend class NegatedMatrix;
		friend class InvertedMatrix;
		friend class RowedMatrix;
		friend class ColedMatrix;
		friend class DiagedMatrix;
		friend class MatedMatrix;
		friend class GetSubMatrix;
		friend class ReturnMatrix;
		friend class LinearEquationSolver;
		friend class GenericMatrix;
		NEW_DELETE(BaseMatrix)
	};


	// ***************************** working classes **************************/

	/// The classes for matrices that can contain data are derived from this.
	class LibGnut_LIBRARY_EXPORT GeneralMatrix : public BaseMatrix         // declarable matrix types
	{
		virtual GeneralMatrix* Image() const;        // copy of matrix
	protected:
		int tag_val;                                 // shows whether can reuse
		int nrows_val, ncols_val;                    // dimensions
		int storage;                                 // total store required
		Real* store;                                 // point to store (0=not set)
		GeneralMatrix();                             // initialise with no store
		GeneralMatrix(ArrayLengthSpecifier);         // constructor getting store
		void Add(GeneralMatrix*, Real);              // sum of GM and Real
		void Add(Real);                              // add Real to this
		void NegAdd(GeneralMatrix*, Real);           // Real - GM
		void NegAdd(Real);                           // this = this - Real
		void Multiply(GeneralMatrix*, Real);         // product of GM and Real
		void Multiply(Real);                         // multiply this by Real
		void Negate(GeneralMatrix*);                 // change sign
		void Negate();                               // change sign
		void ReverseElements();                      // internal reverse of elements
		void ReverseElements(GeneralMatrix*);        // reverse order of elements
		void operator=(Real);                        // set matrix to constant
		Real* GetStore();                            // get store or copy
		GeneralMatrix* BorrowStore(GeneralMatrix*, MatrixType);
		// temporarily access store
		void GetMatrix(const GeneralMatrix*);        // used by = and initialise
		void Eq(const BaseMatrix&, MatrixType);      // used by =
		void Eq(const GeneralMatrix&);               // version with no conversion
		void Eq(const BaseMatrix&, MatrixType, bool);// used by <<
		void Eq2(const BaseMatrix&, MatrixType);     // cut down version of Eq
		int search(const BaseMatrix*) const;
		virtual GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
		void CheckConversion(const BaseMatrix&);     // check conversion OK
		void resize(int, int, int);                  // change dimensions
		virtual short SimpleAddOK(const GeneralMatrix*) { return 0; }
		// see bandmat.cpp for explanation
		virtual void MiniCleanUp()
		{
			store = 0; storage = 0; nrows_val = 0; ncols_val = 0; tag_val = -1;
		}
		// CleanUp when the data array has already been deleted
		void PlusEqual(const GeneralMatrix& gm);
		void MinusEqual(const GeneralMatrix& gm);
		void PlusEqual(Real f);
		void MinusEqual(Real f);
		void swap(GeneralMatrix& gm);                // swap values
	public:
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		virtual MatrixType type() const = 0;         // type of a matrix
		MatrixType Type() const { return type(); }
		int Nrows() const { return nrows_val; }      // get dimensions
		int Ncols() const { return ncols_val; }
		int Storage() const { return storage; }
		Real* Store() const { return store; }
		// updated names
		int nrows() const { return nrows_val; }      // get dimensions
		int ncols() const { return ncols_val; }
		int size() const { return storage; }
		Real* data() { return store; }
		const Real* data() const { return store; }
		const Real* const_data() const { return store; }
		virtual ~GeneralMatrix();                    // delete store if set
		void tDelete();                              // delete if tag_val permits
		bool reuse();                                // true if tag_val allows reuse
		void protect() { tag_val = -1; }               // cannot delete or reuse
		void Protect() { tag_val = -1; }               // cannot delete or reuse
		int tag() const { return tag_val; }
		int Tag() const { return tag_val; }
		bool is_zero() const;                        // test matrix has all zeros
		bool IsZero() const { return is_zero(); }    // test matrix has all zeros
		void Release() { tag_val = 1; }                // del store after next use
		void Release(int t) { tag_val = t; }           // del store after t accesses
		void ReleaseAndDelete() { tag_val = 0; }       // delete matrix after use
		void release() { tag_val = 1; }                // del store after next use
		void release(int t) { tag_val = t; }           // del store after t accesses
		void release_and_delete() { tag_val = 0; }     // delete matrix after use
		void operator<<(const double*);              // assignment from an array
		void operator<<(const float*);               // assignment from an array
		void operator<<(const int*);                 // assignment from an array
		void operator<<(const BaseMatrix& X)
		{
			Eq(X, this->type(), true);
		}              // = without checking type
		void inject(const GeneralMatrix&);           // copy stored els only
		void Inject(const GeneralMatrix& GM_A) { inject(GM_A); }
		void operator+=(const BaseMatrix&);
		void operator-=(const BaseMatrix&);
		void operator*=(const BaseMatrix&);
		void operator|=(const BaseMatrix&);
		void operator&=(const BaseMatrix&);
		void operator+=(Real);
		void operator-=(Real r) { operator+=(-r); }
		void operator*=(Real);
		void operator/=(Real r) { operator*=(1.0 / r); }
		virtual GeneralMatrix* MakeSolver();         // for solving
		virtual void Solver(MatrixColX&, const MatrixColX&) {}
		virtual void GetRow(MatrixRowCol&) = 0;      // Get matrix row
		virtual void RestoreRow(MatrixRowCol&) {}    // Restore matrix row
		virtual void NextRow(MatrixRowCol&);         // Go to next row
		virtual void GetCol(MatrixRowCol&) = 0;      // Get matrix col
		virtual void GetCol(MatrixColX&) = 0;        // Get matrix col
		virtual void RestoreCol(MatrixRowCol&) {}    // Restore matrix col
		virtual void RestoreCol(MatrixColX&) {}      // Restore matrix col
		virtual void NextCol(MatrixRowCol&);         // Go to next col
		virtual void NextCol(MatrixColX&);           // Go to next col
		Real sum_square() const;
		Real sum_absolute_value() const;
		Real sum() const;
		Real maximum_absolute_value1(int& i) const;
		Real minimum_absolute_value1(int& i) const;
		Real maximum1(int& i) const;
		Real minimum1(int& i) const;
		Real maximum_absolute_value() const;
		Real maximum_absolute_value2(int& i, int& j) const;
		Real minimum_absolute_value() const;
		Real minimum_absolute_value2(int& i, int& j) const;
		Real maximum() const;
		Real maximum2(int& i, int& j) const;
		Real minimum() const;
		Real minimum2(int& i, int& j) const;
		LogAndSign log_determinant() const;
		virtual bool IsEqual(const GeneralMatrix&) const;
		// same type, same values
		void CheckStore() const;                     // check store is non-zero
		virtual void SetParameters(const GeneralMatrix*) {}
		// set parameters in GetMatrix
		operator ReturnMatrix() const;               // for building a ReturnMatrix
		ReturnMatrix for_return() const;
		ReturnMatrix ForReturn() const;
		//virtual bool SameStorageType(const GeneralMatrix& A) const;
		//virtual void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
		//virtual void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
		virtual void resize(const GeneralMatrix& A);
		virtual void ReSize(const GeneralMatrix& A) { resize(A); }
		MatrixInput operator<<(double);                // for loading a list
		MatrixInput operator<<(float);                // for loading a list
		MatrixInput operator<<(int f);
		//   ReturnMatrix Reverse() const;                // reverse order of elements
		void cleanup();                              // to clear store

		friend class Matrix;
		friend class SquareMatrix;
		friend class nricMatrix;
		friend class SymmetricMatrix;
		friend class UpperTriangularMatrix;
		friend class LowerTriangularMatrix;
		friend class DiagonalMatrix;
		friend class CroutMatrix;
		friend class RowVector;
		friend class ColumnVector;
		friend class BandMatrix;
		friend class LowerBandMatrix;
		friend class UpperBandMatrix;
		friend class SymmetricBandMatrix;
		friend class BaseMatrix;
		friend class AddedMatrix;
		friend class MultipliedMatrix;
		friend class SubtractedMatrix;
		friend class SPMatrix;
		friend class KPMatrix;
		friend class ConcatenatedMatrix;
		friend class StackedMatrix;
		friend class SolvedMatrix;
		friend class ShiftedMatrix;
		friend class NegShiftedMatrix;
		friend class ScaledMatrix;
		friend class TransposedMatrix;
		friend class ReversedMatrix;
		friend class NegatedMatrix;
		friend class InvertedMatrix;
		friend class RowedMatrix;
		friend class ColedMatrix;
		friend class DiagedMatrix;
		friend class MatedMatrix;
		friend class GetSubMatrix;
		friend class ReturnMatrix;
		friend class LinearEquationSolver;
		friend class GenericMatrix;
		NEW_DELETE(GeneralMatrix)
	};


	/// The usual rectangular matrix.
	class LibGnut_LIBRARY_EXPORT Matrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		Matrix() {}
		~Matrix() {}
		Matrix(int, int);                            // standard declaration
		Matrix(const BaseMatrix&);                   // evaluate BaseMatrix
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const Matrix& m) { Eq(m); }
		MatrixType type() const;
		Real& operator()(int, int);                  // access element
		Real& element(int, int);                     // access element
		Real operator()(int, int) const;             // access element
		Real element(int, int) const;                // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + (int64_t)m * ncols_val; }
		const Real* operator[](int m) const { return store + (int64_t)m * ncols_val; }
		// following for Numerical Recipes in C++
		Matrix(Real, int, int);
		Matrix(const Real*, int, int);
#endif
		Matrix(const Matrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
		GeneralMatrix* MakeSolver();
		Real trace() const;
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void RestoreCol(MatrixRowCol&);
		void RestoreCol(MatrixColX&);
		void NextRow(MatrixRowCol&);
		void NextCol(MatrixRowCol&);
		void NextCol(MatrixColX&);
		virtual void resize(int, int);           // change dimensions
		   // virtual so we will catch it being used in a vector called as a matrix
		virtual void resize_keep(int, int);
		virtual void ReSize(int m, int n) { resize(m, n); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		Real maximum_absolute_value2(int& i, int& j) const;
		Real minimum_absolute_value2(int& i, int& j) const;
		Real maximum2(int& i, int& j) const;
		Real minimum2(int& i, int& j) const;
		void operator+=(const Matrix& M) { PlusEqual(M); }
		void operator-=(const Matrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::Add(f); }
		void operator-=(Real f) { GeneralMatrix::Add(-f); }
		void swap(Matrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
		LibGnut_LIBRARY_EXPORT friend Real dotproduct(const Matrix& A, const Matrix& B);
		NEW_DELETE(Matrix)
	};

	/// Square matrix.
	class LibGnut_LIBRARY_EXPORT SquareMatrix : public Matrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		SquareMatrix() {}
		~SquareMatrix() {}
		SquareMatrix(ArrayLengthSpecifier);          // standard declaration
		SquareMatrix(const BaseMatrix&);             // evaluate BaseMatrix
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const SquareMatrix& m) { Eq(m); }
		void operator=(const Matrix& m);
		MatrixType type() const;
		SquareMatrix(const SquareMatrix& gm) : Matrix() { GetMatrix(&gm); }
		SquareMatrix(const Matrix& gm);
		void resize(int);                            // change dimensions
		void ReSize(int m) { resize(m); }
		void resize_keep(int);
		void resize_keep(int, int);
		void resize(int, int);                        // change dimensions
		void ReSize(int m, int n) { resize(m, n); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		void operator+=(const Matrix& M) { PlusEqual(M); }
		void operator-=(const Matrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::Add(f); }
		void operator-=(Real f) { GeneralMatrix::Add(-f); }
		void swap(SquareMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
		NEW_DELETE(SquareMatrix)
	};

	/// Rectangular matrix for use with Numerical Recipes in C.
	class nricMatrix : public Matrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
		Real** row_pointer;                          // points to rows
		void MakeRowPointer();                       // build rowpointer
		void DeleteRowPointer();
	public:
		nricMatrix() {}
		nricMatrix(int m, int n)                     // standard declaration
			: Matrix(m, n) {
			MakeRowPointer();
		}
		nricMatrix(const BaseMatrix& bm)             // evaluate BaseMatrix
			: Matrix(bm) {
			MakeRowPointer();
		}
		void operator=(const BaseMatrix& bm)
		{
			DeleteRowPointer(); Matrix::operator=(bm); MakeRowPointer();
		}
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const nricMatrix& m)
		{
			DeleteRowPointer(); Eq(m); MakeRowPointer();
		}
		void operator<<(const BaseMatrix& X)
		{
			DeleteRowPointer(); Eq(X, this->type(), true); MakeRowPointer();
		}
		nricMatrix(const nricMatrix& gm) : Matrix()
		{
			GetMatrix(&gm); MakeRowPointer();
		}
		void resize(int m, int n)               // change dimensions
		{
			DeleteRowPointer(); Matrix::resize(m, n); MakeRowPointer();
		}
		void resize_keep(int m, int n)               // change dimensions
		{
			DeleteRowPointer(); Matrix::resize_keep(m, n); MakeRowPointer();
		}
		void ReSize(int m, int n)               // change dimensions
		{
			DeleteRowPointer(); Matrix::resize(m, n); MakeRowPointer();
		}
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		~nricMatrix() { DeleteRowPointer(); }
		Real** nric() const { CheckStore(); return row_pointer - 1; }
		void cleanup();                                // to clear store
		void MiniCleanUp();
		void operator+=(const Matrix& M) { PlusEqual(M); }
		void operator-=(const Matrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::Add(f); }
		void operator-=(Real f) { GeneralMatrix::Add(-f); }
		void swap(nricMatrix& gm);
		NEW_DELETE(nricMatrix)
	};

	/// Symmetric matrix.
	class LibGnut_LIBRARY_EXPORT SymmetricMatrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		SymmetricMatrix() {}
		~SymmetricMatrix() {}
		SymmetricMatrix(ArrayLengthSpecifier);
		SymmetricMatrix(const BaseMatrix&);
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const SymmetricMatrix& m) { Eq(m); }
		Real& operator()(int, int);                  // access element
		Real& element(int, int);                     // access element
		Real operator()(int, int) const;             // access element
		Real element(int, int) const;                // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + (m * (m + 1)) / 2; }
		const Real* operator[](int m) const { return store + (m * (m + 1)) / 2; }
#endif
		MatrixType type() const;
		SymmetricMatrix(const SymmetricMatrix& gm)
			: GeneralMatrix() {
			GetMatrix(&gm);
		}
		Real sum_square() const;
		Real sum_absolute_value() const;
		Real sum() const;
		Real trace() const;
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void RestoreCol(MatrixRowCol&) {}
		void RestoreCol(MatrixColX&);
		GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
		void resize(int);                           // change dimensions
		void ReSize(int m) { resize(m); }
		void resize_keep(int);
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		void operator+=(const SymmetricMatrix& M) { PlusEqual(M); }
		void operator-=(const SymmetricMatrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::Add(f); }
		void operator-=(Real f) { GeneralMatrix::Add(-f); }
		void swap(SymmetricMatrix& gm) { GeneralMatrix::swap((GeneralMatrix&)gm); }
		NEW_DELETE(SymmetricMatrix)
	};

	/// Upper triangular matrix.
	class LibGnut_LIBRARY_EXPORT UpperTriangularMatrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		UpperTriangularMatrix() {}
		~UpperTriangularMatrix() {}
		UpperTriangularMatrix(ArrayLengthSpecifier);
		void operator=(const BaseMatrix&);
		void operator=(const UpperTriangularMatrix& m) { Eq(m); }
		UpperTriangularMatrix(const BaseMatrix&);
		UpperTriangularMatrix(const UpperTriangularMatrix& gm)
			: GeneralMatrix() {
			GetMatrix(&gm);
		}
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		Real& operator()(int, int);                  // access element
		Real& element(int, int);                     // access element
		Real operator()(int, int) const;             // access element
		Real element(int, int) const;                // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + (int64_t)m * ncols_val - ((int64_t)m * ((int64_t)m + 1)) / 2; }
		const Real* operator[](int m) const
		{
			return store + (int64_t)m * ncols_val - ((int64_t)m * ((int64_t)m + 1)) / 2;
		}
#endif
		MatrixType type() const;
		GeneralMatrix* MakeSolver() { return this; } // for solving
		void Solver(MatrixColX&, const MatrixColX&);
		LogAndSign log_determinant() const;
		Real trace() const;
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void RestoreCol(MatrixRowCol&);
		void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
		void NextRow(MatrixRowCol&);
		void resize(int);                       // change dimensions
		void ReSize(int m) { resize(m); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		void resize_keep(int);
		MatrixBandWidth bandwidth() const;
		void operator+=(const UpperTriangularMatrix& M) { PlusEqual(M); }
		void operator-=(const UpperTriangularMatrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::operator+=(f); }
		void operator-=(Real f) { GeneralMatrix::operator-=(f); }
		void swap(UpperTriangularMatrix& gm)
		{
			GeneralMatrix::swap((GeneralMatrix&)gm);
		}
		NEW_DELETE(UpperTriangularMatrix)
	};

	/// Lower triangular matrix.
	class LibGnut_LIBRARY_EXPORT LowerTriangularMatrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		LowerTriangularMatrix() {}
		~LowerTriangularMatrix() {}
		LowerTriangularMatrix(ArrayLengthSpecifier);
		LowerTriangularMatrix(const LowerTriangularMatrix& gm)
			: GeneralMatrix() {
			GetMatrix(&gm);
		}
		LowerTriangularMatrix(const BaseMatrix& M);
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const LowerTriangularMatrix& m) { Eq(m); }
		Real& operator()(int, int);                  // access element
		Real& element(int, int);                     // access element
		Real operator()(int, int) const;             // access element
		Real element(int, int) const;                // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + (m * (m + 1)) / 2; }
		const Real* operator[](int m) const { return store + (m * (m + 1)) / 2; }
#endif
		MatrixType type() const;
		GeneralMatrix* MakeSolver() { return this; } // for solving
		void Solver(MatrixColX&, const MatrixColX&);
		LogAndSign log_determinant() const;
		Real trace() const;
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void RestoreCol(MatrixRowCol&);
		void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
		void NextRow(MatrixRowCol&);
		void resize(int);                       // change dimensions
		void ReSize(int m) { resize(m); }
		void resize_keep(int);
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		MatrixBandWidth bandwidth() const;
		void operator+=(const LowerTriangularMatrix& M) { PlusEqual(M); }
		void operator-=(const LowerTriangularMatrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::operator+=(f); }
		void operator-=(Real f) { GeneralMatrix::operator-=(f); }
		void swap(LowerTriangularMatrix& gm)
		{
			GeneralMatrix::swap((GeneralMatrix&)gm);
		}
		NEW_DELETE(LowerTriangularMatrix)
	};

	/// Diagonal matrix.
	class LibGnut_LIBRARY_EXPORT DiagonalMatrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		DiagonalMatrix() {}
		~DiagonalMatrix() {}
		DiagonalMatrix(ArrayLengthSpecifier);
		DiagonalMatrix(const BaseMatrix&);
		DiagonalMatrix(const DiagonalMatrix& gm)
			: GeneralMatrix() {
			GetMatrix(&gm);
		}
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const DiagonalMatrix& m) { Eq(m); }
		Real& operator()(int, int);                  // access element
		Real& operator()(int);                       // access element
		Real operator()(int, int) const;             // access element
		Real operator()(int) const;
		Real& element(int, int);                     // access element
		Real& element(int);                          // access element
		Real element(int, int) const;                // access element
		Real element(int) const;                     // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real& operator[](int m) { return store[m]; }
		const Real& operator[](int m) const { return store[m]; }
#endif
		MatrixType type() const;

		LogAndSign log_determinant() const;
		Real trace() const;
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void NextRow(MatrixRowCol&);
		void NextCol(MatrixRowCol&);
		void NextCol(MatrixColX&);
		GeneralMatrix* MakeSolver() { return this; } // for solving
		void Solver(MatrixColX&, const MatrixColX&);
		GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
		void resize(int);                       // change dimensions
		void ReSize(int m) { resize(m); }
		void resize_keep(int);
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		Real* nric() const
		{
			CheckStore(); return store - 1;
		}         // for use by NRIC
		MatrixBandWidth bandwidth() const;
		//   ReturnMatrix Reverse() const;                // reverse order of elements
		void operator+=(const DiagonalMatrix& M) { PlusEqual(M); }
		void operator-=(const DiagonalMatrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::operator+=(f); }
		void operator-=(Real f) { GeneralMatrix::operator-=(f); }
		void swap(DiagonalMatrix& gm)
		{
			GeneralMatrix::swap((GeneralMatrix&)gm);
		}
		NEW_DELETE(DiagonalMatrix)
	};

	/// Row vector.
	class LibGnut_LIBRARY_EXPORT RowVector : public Matrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		RowVector() { nrows_val = 1; }
		~RowVector() {}
		RowVector(ArrayLengthSpecifier n) : Matrix(1, n.Value()) {}
		RowVector(const BaseMatrix&);
		RowVector(const RowVector& gm) : Matrix() { GetMatrix(&gm); }
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const RowVector& m) { Eq(m); }
		Real& operator()(int);                       // access element
		Real& element(int);                          // access element
		Real operator()(int) const;                  // access element
		Real element(int) const;                     // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real& operator[](int m) { return store[m]; }
		const Real& operator[](int m) const { return store[m]; }
		// following for Numerical Recipes in C++
		RowVector(Real a, int n) : Matrix(a, 1, n) {}
		RowVector(const Real* a, int n) : Matrix(a, 1, n) {}
#endif
		MatrixType type() const;
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void NextCol(MatrixRowCol&);
		void NextCol(MatrixColX&);
		void RestoreCol(MatrixRowCol&) {}
		void RestoreCol(MatrixColX& c);
		GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
		void resize(int);                       // change dimensions
		void ReSize(int m) { resize(m); }
		void resize_keep(int);
		void resize_keep(int, int);
		void resize(int, int);                   // in case access is matrix
		void ReSize(int m, int n) { resize(m, n); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		Real* nric() const
		{
			CheckStore(); return store - 1;
		}         // for use by NRIC
		void cleanup();                              // to clear store
		void MiniCleanUp()
		{
			store = 0; storage = 0; nrows_val = 1; ncols_val = 0; tag_val = -1;
		}
		// friend ReturnMatrix GetMatrixRow(Matrix& A, int row);
		void operator+=(const Matrix& M) { PlusEqual(M); }
		void operator-=(const Matrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::Add(f); }
		void operator-=(Real f) { GeneralMatrix::Add(-f); }
		void swap(RowVector& gm)
		{
			GeneralMatrix::swap((GeneralMatrix&)gm);
		}
		NEW_DELETE(RowVector)
	};

	/// Column vector.
	class LibGnut_LIBRARY_EXPORT ColumnVector : public Matrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		ColumnVector() { ncols_val = 1; }
		~ColumnVector() {}
		ColumnVector(ArrayLengthSpecifier n) : Matrix(n.Value(), 1) {}
		ColumnVector(const BaseMatrix&);
		ColumnVector(const ColumnVector& gm) : Matrix() { GetMatrix(&gm); }
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const ColumnVector& m) { Eq(m); }
		Real& operator()(int);                       // access element
		Real& element(int);                          // access element
		Real operator()(int) const;                  // access element
		Real element(int) const;                     // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real& operator[](int m) { return store[m]; }
		const Real& operator[](int m) const { return store[m]; }
		// following for Numerical Recipes in C++
		ColumnVector(Real a, int m) : Matrix(a, m, 1) {}
		ColumnVector(const Real* a, int m) : Matrix(a, m, 1) {}
#endif
		MatrixType type() const;
		GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
		void resize(int);                       // change dimensions
		void ReSize(int m) { resize(m); }
		void resize_keep(int);
		void resize_keep(int, int);
		void resize(int, int);                   // in case access is matrix
		void ReSize(int m, int n) { resize(m, n); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		Real* nric() const
		{
			CheckStore(); return store - 1;
		}         // for use by NRIC
		void cleanup();                              // to clear store
		void MiniCleanUp()
		{
			store = 0; storage = 0; nrows_val = 0; ncols_val = 1; tag_val = -1;
		}
		//   ReturnMatrix Reverse() const;                // reverse order of elements
		void operator+=(const Matrix& M) { PlusEqual(M); }
		void operator-=(const Matrix& M) { MinusEqual(M); }
		void operator+=(Real f) { GeneralMatrix::Add(f); }
		void operator-=(Real f) { GeneralMatrix::Add(-f); }
		void swap(ColumnVector& gm)
		{
			GeneralMatrix::swap((GeneralMatrix&)gm);
		}
		NEW_DELETE(ColumnVector)
	};

	/// LU matrix.
	/// A square matrix decomposed into upper and lower triangular
	/// in preparation for inverting or solving equations.
	class LibGnut_LIBRARY_EXPORT CroutMatrix : public GeneralMatrix
	{
		int* indx;
		bool d;                              // number of row swaps = even or odd
		bool sing;
		void ludcmp();
		void get_aux(CroutMatrix&);                  // for copying indx[] etc
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		CroutMatrix(const BaseMatrix&);
		CroutMatrix() : indx(0), d(true), sing(true) {}
		CroutMatrix(const CroutMatrix&);
		void operator=(const CroutMatrix&);
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixType type() const;
		void lubksb(Real*, int = 0);
		~CroutMatrix();
		GeneralMatrix* MakeSolver() { return this; } // for solving
		LogAndSign log_determinant() const;
		void Solver(MatrixColX&, const MatrixColX&);
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
		void cleanup();                                // to clear store
		void MiniCleanUp();
		bool IsEqual(const GeneralMatrix&) const;
		bool is_singular() const { return sing; }
		bool IsSingular() const { return sing; }
		const int* const_data_indx() const { return indx; }
		bool even_exchanges() const { return d; }
		void swap(CroutMatrix& gm);
		NEW_DELETE(CroutMatrix)
	};

	// ***************************** band matrices ***************************/

	/// Band matrix.
	class LibGnut_LIBRARY_EXPORT BandMatrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	protected:
		void CornerClear() const;                    // set unused elements to zero
		short SimpleAddOK(const GeneralMatrix* gm);
	public:
		int lower_val, upper_val;                            // band widths
		BandMatrix() { lower_val = 0; upper_val = 0; CornerClear(); }
		~BandMatrix() {}
		BandMatrix(int n, int lb, int ub) { resize(n, lb, ub); CornerClear(); }
		// standard declaration
		BandMatrix(const BaseMatrix&);               // evaluate BaseMatrix
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const BandMatrix& m) { Eq(m); }
		MatrixType type() const;
		Real& operator()(int, int);                  // access element
		Real& element(int, int);                     // access element
		Real operator()(int, int) const;             // access element
		Real element(int, int) const;                // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + ((int64_t)upper_val + lower_val) * (int64_t)m + lower_val; }
		const Real* operator[](int m) const
		{
			return store + ((int64_t)upper_val + lower_val) * (int64_t)m + lower_val;
		}
#endif
		BandMatrix(const BandMatrix& gm) : GeneralMatrix() { GetMatrix(&gm); }
		LogAndSign log_determinant() const;
		GeneralMatrix* MakeSolver();
		Real trace() const;
		Real sum_square() const
		{
			CornerClear(); return GeneralMatrix::sum_square();
		}
		Real sum_absolute_value() const
		{
			CornerClear(); return GeneralMatrix::sum_absolute_value();
		}
		Real sum() const
		{
			CornerClear(); return GeneralMatrix::sum();
		}
		Real maximum_absolute_value() const
		{
			CornerClear(); return GeneralMatrix::maximum_absolute_value();
		}
		Real minimum_absolute_value() const
		{
			int i, j; return GeneralMatrix::minimum_absolute_value2(i, j);
		}
		Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
		Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void RestoreCol(MatrixRowCol&);
		void RestoreCol(MatrixColX& c) { RestoreCol((MatrixRowCol&)c); }
		void NextRow(MatrixRowCol&);
		virtual void resize(int, int, int);             // change dimensions
		virtual void ReSize(int m, int n, int b) { resize(m, n, b); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		//bool SameStorageType(const GeneralMatrix& A) const;
		//void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
		//void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
		MatrixBandWidth bandwidth() const;
		void SetParameters(const GeneralMatrix*);
		MatrixInput operator<<(double);                // will give error
		MatrixInput operator<<(float);                // will give error
		MatrixInput operator<<(int f);
		void operator<<(const double* r);              // will give error
		void operator<<(const float* r);              // will give error
		void operator<<(const int* r);               // will give error
		   // the next is included because Zortech and Borland
		   // cannot find the copy in GeneralMatrix
		void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
		void swap(BandMatrix& gm);
		NEW_DELETE(BandMatrix)
	};

	/// Upper triangular band matrix.
	class LibGnut_LIBRARY_EXPORT UpperBandMatrix : public BandMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		UpperBandMatrix() {}
		~UpperBandMatrix() {}
		UpperBandMatrix(int n, int ubw)              // standard declaration
			: BandMatrix(n, 0, ubw) {}
		UpperBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const UpperBandMatrix& m) { Eq(m); }
		MatrixType type() const;
		UpperBandMatrix(const UpperBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
		GeneralMatrix* MakeSolver() { return this; }
		void Solver(MatrixColX&, const MatrixColX&);
		LogAndSign log_determinant() const;
		void resize(int, int, int);             // change dimensions
		void ReSize(int m, int n, int b) { resize(m, n, b); }
		void resize(int n, int ubw)              // change dimensions
		{
			BandMatrix::resize(n, 0, ubw);
		}
		void ReSize(int n, int ubw)              // change dimensions
		{
			BandMatrix::resize(n, 0, ubw);
		}
		void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
		void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
		Real& operator()(int, int);
		Real operator()(int, int) const;
		Real& element(int, int);
		Real element(int, int) const;
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + upper_val * (int64_t)m; }
		const Real* operator[](int m) const { return store + upper_val * (int64_t)m; }
#endif
		void swap(UpperBandMatrix& gm)
		{
			BandMatrix::swap((BandMatrix&)gm);
		}
		NEW_DELETE(UpperBandMatrix)
	};

	/// Lower triangular band matrix.
	class LibGnut_LIBRARY_EXPORT LowerBandMatrix : public BandMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		LowerBandMatrix() {}
		~LowerBandMatrix() {}
		LowerBandMatrix(int n, int lbw)              // standard declaration
			: BandMatrix(n, lbw, 0) {}
		LowerBandMatrix(const BaseMatrix&);          // evaluate BaseMatrix
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const LowerBandMatrix& m) { Eq(m); }
		MatrixType type() const;
		LowerBandMatrix(const LowerBandMatrix& gm) : BandMatrix() { GetMatrix(&gm); }
		GeneralMatrix* MakeSolver() { return this; }
		void Solver(MatrixColX&, const MatrixColX&);
		LogAndSign log_determinant() const;
		void resize(int, int, int);             // change dimensions
		void ReSize(int m, int n, int b) { resize(m, n, b); }
		void resize(int n, int lbw)             // change dimensions
		{
			BandMatrix::resize(n, lbw, 0);
		}
		void ReSize(int n, int lbw)             // change dimensions
		{
			BandMatrix::resize(n, lbw, 0);
		}
		void resize(const GeneralMatrix& A) { BandMatrix::resize(A); }
		void ReSize(const GeneralMatrix& A) { BandMatrix::resize(A); }
		Real& operator()(int, int);
		Real operator()(int, int) const;
		Real& element(int, int);
		Real element(int, int) const;
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + lower_val * ((int64_t)m + 1); }
		const Real* operator[](int m) const { return store + lower_val * ((int64_t)m + 1); }
#endif
		void swap(LowerBandMatrix& gm)
		{
			BandMatrix::swap((BandMatrix&)gm);
		}
		NEW_DELETE(LowerBandMatrix)
	};

	/// Symmetric band matrix.
	class LibGnut_LIBRARY_EXPORT SymmetricBandMatrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;                // copy of matrix
		void CornerClear() const;                    // set unused elements to zero
		short SimpleAddOK(const GeneralMatrix* gm);
	public:
		int lower_val;                                   // lower band width
		SymmetricBandMatrix() { lower_val = 0; CornerClear(); }
		~SymmetricBandMatrix() {}
		SymmetricBandMatrix(int n, int lb) { resize(n, lb); CornerClear(); }
		SymmetricBandMatrix(const BaseMatrix&);
		void operator=(const BaseMatrix&);
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		void operator=(const SymmetricBandMatrix& m) { Eq(m); }
		Real& operator()(int, int);                  // access element
		Real& element(int, int);                     // access element
		Real operator()(int, int) const;             // access element
		Real element(int, int) const;                // access element
#ifdef SETUP_C_SUBSCRIPTS
		Real* operator[](int m) { return store + lower_val * ((int64_t)m + 1); }
		const Real* operator[](int m) const { return store + lower_val * ((int64_t)m + 1); }
#endif
		MatrixType type() const;
		SymmetricBandMatrix(const SymmetricBandMatrix& gm)
			: GeneralMatrix() {
			GetMatrix(&gm);
		}
		GeneralMatrix* MakeSolver();
		Real sum_square() const;
		Real sum_absolute_value() const;
		Real sum() const;
		Real maximum_absolute_value() const
		{
			CornerClear(); return GeneralMatrix::maximum_absolute_value();
		}
		Real minimum_absolute_value() const
		{
			int i, j; return GeneralMatrix::minimum_absolute_value2(i, j);
		}
		Real maximum() const { int i, j; return GeneralMatrix::maximum2(i, j); }
		Real minimum() const { int i, j; return GeneralMatrix::minimum2(i, j); }
		Real trace() const;
		LogAndSign log_determinant() const;
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void RestoreCol(MatrixRowCol&) {}
		void RestoreCol(MatrixColX&);
		GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
		void resize(int, int);                       // change dimensions
		void ReSize(int m, int b) { resize(m, b); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		//bool SameStorageType(const GeneralMatrix& A) const;
		//void ReSizeForAdd(const GeneralMatrix& A, const GeneralMatrix& B);
		//void ReSizeForSP(const GeneralMatrix& A, const GeneralMatrix& B);
		MatrixBandWidth bandwidth() const;
		void SetParameters(const GeneralMatrix*);
		void operator<<(const double* r);              // will give error
		void operator<<(const float* r);              // will give error
		void operator<<(const int* r);               // will give error
		void operator<<(const BaseMatrix& X) { GeneralMatrix::operator<<(X); }
		void swap(SymmetricBandMatrix& gm);
		NEW_DELETE(SymmetricBandMatrix)
	};

	/// LU decomposition of a band matrix.
	class LibGnut_LIBRARY_EXPORT BandLUMatrix : public GeneralMatrix
	{
		int* indx;
		bool d;
		bool sing;                                   // true if singular
		Real* store2;
		int storage2;
		int m1, m2;                                   // lower and upper
		void ludcmp();
		void get_aux(BandLUMatrix&);                 // for copying indx[] etc
		GeneralMatrix* Image() const;                // copy of matrix
	public:
		BandLUMatrix(const BaseMatrix&);
		BandLUMatrix()
			: indx(0), d(true), sing(true), store2(0), storage2(0), m1(0), m2(0) {}
		BandLUMatrix(const BandLUMatrix&);
		void operator=(const BandLUMatrix&);
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixType type() const;
		void lubksb(Real*, int = 0);
		~BandLUMatrix();
		GeneralMatrix* MakeSolver() { return this; } // for solving
		LogAndSign log_determinant() const;
		void Solver(MatrixColX&, const MatrixColX&);
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX& c) { GetCol((MatrixRowCol&)c); }
		void cleanup();                                // to clear store
		void MiniCleanUp();
		bool IsEqual(const GeneralMatrix&) const;
		bool is_singular() const { return sing; }
		bool IsSingular() const { return sing; }
		const Real* const_data2() const { return store2; }
		int size2() const { return storage2; }
		const int* const_data_indx() const { return indx; }
		bool even_exchanges() const { return d; }
		MatrixBandWidth bandwidth() const;
		void swap(BandLUMatrix& gm);
		NEW_DELETE(BandLUMatrix)
	};

	// ************************** special matrices ****************************

	/// Identity matrix.
	class LibGnut_LIBRARY_EXPORT IdentityMatrix : public GeneralMatrix
	{
		GeneralMatrix* Image() const;          // copy of matrix
	public:
		IdentityMatrix() {}
		~IdentityMatrix() {}
		IdentityMatrix(ArrayLengthSpecifier n) : GeneralMatrix(1)
		{
			nrows_val = ncols_val = n.Value(); *store = 1;
		}
		IdentityMatrix(const IdentityMatrix& gm)
			: GeneralMatrix() {
			GetMatrix(&gm);
		}
		IdentityMatrix(const BaseMatrix&);
		void operator=(const BaseMatrix&);
		void operator=(const IdentityMatrix& m) { Eq(m); }
		void operator=(Real f) { GeneralMatrix::operator=(f); }
		MatrixType type() const;

		LogAndSign log_determinant() const;
		Real trace() const;
		Real sum_square() const;
		Real sum_absolute_value() const;
		Real sum() const { return trace(); }
		void GetRow(MatrixRowCol&);
		void GetCol(MatrixRowCol&);
		void GetCol(MatrixColX&);
		void NextRow(MatrixRowCol&);
		void NextCol(MatrixRowCol&);
		void NextCol(MatrixColX&);
		GeneralMatrix* MakeSolver() { return this; } // for solving
		void Solver(MatrixColX&, const MatrixColX&);
		GeneralMatrix* Transpose(TransposedMatrix*, MatrixType);
		void resize(int n);
		void ReSize(int n) { resize(n); }
		void resize(const GeneralMatrix& A);
		void ReSize(const GeneralMatrix& A) { resize(A); }
		MatrixBandWidth bandwidth() const;
		//   ReturnMatrix Reverse() const;                // reverse order of elements
		void swap(IdentityMatrix& gm)
		{
			GeneralMatrix::swap((GeneralMatrix&)gm);
		}
		NEW_DELETE(IdentityMatrix)
	};




	// ************************** GenericMatrix class ************************/

	/// A matrix which can be of any GeneralMatrix type.
	class LibGnut_LIBRARY_EXPORT GenericMatrix : public BaseMatrix
	{
		GeneralMatrix* gm;
		int search(const BaseMatrix* bm) const;
		friend class BaseMatrix;
	public:
		GenericMatrix() : gm(0) {}
		GenericMatrix(const BaseMatrix& bm)
		{
			gm = ((BaseMatrix&)bm).Evaluate(); gm = gm->Image();
		}
		GenericMatrix(const GenericMatrix& bm) : BaseMatrix()
		{
			gm = bm.gm->Image();
		}
		void operator=(const GenericMatrix&);
		void operator=(const BaseMatrix&);
		void operator+=(const BaseMatrix&);
		void operator-=(const BaseMatrix&);
		void operator*=(const BaseMatrix&);
		void operator|=(const BaseMatrix&);
		void operator&=(const BaseMatrix&);
		void operator+=(Real);
		void operator-=(Real r) { operator+=(-r); }
		void operator*=(Real);
		void operator/=(Real r) { operator*=(1.0 / r); }
		~GenericMatrix() { delete gm; }
		void cleanup() { delete gm; gm = 0; }
		void Release() { gm->Release(); }
		void release() { gm->release(); }
		GeneralMatrix* Evaluate(MatrixType = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		void swap(GenericMatrix& x);
		NEW_DELETE(GenericMatrix)
	};

	// *************************** temporary classes *************************/

	/// Product of two matrices.
	/// \internal
	class LibGnut_LIBRARY_EXPORT MultipliedMatrix : public BaseMatrix
	{
	protected:
		// if these union statements cause problems, simply remove them
		// and declare the items individually
		union { const BaseMatrix* bm1; GeneralMatrix* gm1; };
		// pointers to summands
		union { const BaseMatrix* bm2; GeneralMatrix* gm2; };
		MultipliedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: bm1(bm1x), bm2(bm2x) {}
		int search(const BaseMatrix*) const;
		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~MultipliedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(MultipliedMatrix)
	};

	/// Sum of two matrices.
	/// \internal
	class LibGnut_LIBRARY_EXPORT AddedMatrix : public MultipliedMatrix
	{
	protected:
		AddedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: MultipliedMatrix(bm1x, bm2x) {}

		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~AddedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(AddedMatrix)
	};

	/// Schur (elementwise) product of two matrices.
	/// \internal
	class LibGnut_LIBRARY_EXPORT SPMatrix : public AddedMatrix
	{
	protected:
		SPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: AddedMatrix(bm1x, bm2x) {}

		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~SPMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;

		LibGnut_LIBRARY_EXPORT friend SPMatrix SP(const BaseMatrix&, const BaseMatrix&);

		NEW_DELETE(SPMatrix)
	};

	/// Kronecker product of two matrices.
	/// \internal
	class LibGnut_LIBRARY_EXPORT KPMatrix : public MultipliedMatrix
	{
	protected:
		KPMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: MultipliedMatrix(bm1x, bm2x) {}

		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~KPMatrix() {}
		MatrixBandWidth bandwidth() const;
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		LibGnut_LIBRARY_EXPORT friend KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
		NEW_DELETE(KPMatrix)
	};

	/// Two matrices horizontally concatenated.
	/// \internal
	class LibGnut_LIBRARY_EXPORT ConcatenatedMatrix : public MultipliedMatrix
	{
	protected:
		ConcatenatedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: MultipliedMatrix(bm1x, bm2x) {}

		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~ConcatenatedMatrix() {}
		MatrixBandWidth bandwidth() const;
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		NEW_DELETE(ConcatenatedMatrix)
	};

	/// Two matrices vertically concatenated.
	/// \internal
	class LibGnut_LIBRARY_EXPORT StackedMatrix : public ConcatenatedMatrix
	{
	protected:
		StackedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: ConcatenatedMatrix(bm1x, bm2x) {}

		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~StackedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		NEW_DELETE(StackedMatrix)
	};

	/// Inverted matrix times matrix.
	/// \internal
	class LibGnut_LIBRARY_EXPORT SolvedMatrix : public MultipliedMatrix
	{
		SolvedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: MultipliedMatrix(bm1x, bm2x) {}
		friend class BaseMatrix;
		friend class InvertedMatrix;                        // for operator*
	public:
		~SolvedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(SolvedMatrix)
	};

	/// Difference between two matrices.
	/// \internal
	class LibGnut_LIBRARY_EXPORT SubtractedMatrix : public AddedMatrix
	{
		SubtractedMatrix(const BaseMatrix* bm1x, const BaseMatrix* bm2x)
			: AddedMatrix(bm1x, bm2x) {}
		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~SubtractedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		NEW_DELETE(SubtractedMatrix)
	};

	/// Any type of matrix plus Real.
	/// \internal
	class LibGnut_LIBRARY_EXPORT ShiftedMatrix : public BaseMatrix
	{
	protected:
		union { const BaseMatrix* bm; GeneralMatrix* gm; };
		Real f;
		ShiftedMatrix(const BaseMatrix* bmx, Real fx) : bm(bmx), f(fx) {}
		int search(const BaseMatrix*) const;
		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~ShiftedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		LibGnut_LIBRARY_EXPORT friend ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
		NEW_DELETE(ShiftedMatrix)
	};

	/// Real minus matrix.
	/// \internal
	class LibGnut_LIBRARY_EXPORT NegShiftedMatrix : public ShiftedMatrix
	{
	protected:
		NegShiftedMatrix(Real fx, const BaseMatrix* bmx) : ShiftedMatrix(bmx, fx) {}
		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~NegShiftedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		LibGnut_LIBRARY_EXPORT friend NegShiftedMatrix operator-(Real, const BaseMatrix&);
		NEW_DELETE(NegShiftedMatrix)
	};

	/// Any type of matrix times Real.
	/// \internal
	class LibGnut_LIBRARY_EXPORT ScaledMatrix : public ShiftedMatrix
	{
		ScaledMatrix(const BaseMatrix* bmx, Real fx) : ShiftedMatrix(bmx, fx) {}
		friend class BaseMatrix;
		friend class GeneralMatrix;
		friend class GenericMatrix;
	public:
		~ScaledMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		LibGnut_LIBRARY_EXPORT friend ScaledMatrix operator*(Real f, const BaseMatrix& BM);
		NEW_DELETE(ScaledMatrix)
	};

	/// Any type of matrix times -1.
	/// \internal
	class LibGnut_LIBRARY_EXPORT NegatedMatrix : public BaseMatrix
	{
	protected:
		union { const BaseMatrix* bm; GeneralMatrix* gm; };
		NegatedMatrix(const BaseMatrix* bmx) : bm(bmx) {}
		int search(const BaseMatrix*) const;
	private:
		friend class BaseMatrix;
	public:
		~NegatedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(NegatedMatrix)
	};

	/// Transposed matrix.
	/// \internal
	class LibGnut_LIBRARY_EXPORT TransposedMatrix : public NegatedMatrix
	{
		TransposedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
		friend class BaseMatrix;
	public:
		~TransposedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(TransposedMatrix)
	};

	/// Any type of matrix with order of elements reversed.
	/// \internal
	class LibGnut_LIBRARY_EXPORT ReversedMatrix : public NegatedMatrix
	{
		ReversedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
		friend class BaseMatrix;
	public:
		~ReversedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		NEW_DELETE(ReversedMatrix)
	};

	/// Inverse of matrix.
	/// \internal
	class LibGnut_LIBRARY_EXPORT InvertedMatrix : public NegatedMatrix
	{
		InvertedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
	public:
		~InvertedMatrix() {}
		SolvedMatrix operator*(const BaseMatrix&) const;       // inverse(A) * B
		ScaledMatrix operator*(Real t) const { return BaseMatrix::operator*(t); }
		friend class BaseMatrix;
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(InvertedMatrix)
	};

	/// Any type of matrix interpreted as a RowVector.
	/// \internal
	class LibGnut_LIBRARY_EXPORT RowedMatrix : public NegatedMatrix
	{
		RowedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
		friend class BaseMatrix;
	public:
		~RowedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(RowedMatrix)
	};

	/// Any type of matrix interpreted as a ColumnVector.
	/// \internal
	class LibGnut_LIBRARY_EXPORT ColedMatrix : public NegatedMatrix
	{
		ColedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
		friend class BaseMatrix;
	public:
		~ColedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(ColedMatrix)
	};

	/// Any type of matrix interpreted as a DiagonalMatrix.
	/// \internal
	class LibGnut_LIBRARY_EXPORT DiagedMatrix : public NegatedMatrix
	{
		DiagedMatrix(const BaseMatrix* bmx) : NegatedMatrix(bmx) {}
		friend class BaseMatrix;
	public:
		~DiagedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(DiagedMatrix)
	};

	/// Any type of matrix interpreted as a (rectangular) Matrix.
	/// \internal
	class LibGnut_LIBRARY_EXPORT MatedMatrix : public NegatedMatrix
	{
		int nr, nc;
		MatedMatrix(const BaseMatrix* bmx, int nrx, int ncx)
			: NegatedMatrix(bmx), nr(nrx), nc(ncx) {}
		friend class BaseMatrix;
	public:
		~MatedMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(MatedMatrix)
	};

	/// A matrix in an "envelope' for return from a function.
	/// \internal
	class LibGnut_LIBRARY_EXPORT ReturnMatrix : public BaseMatrix
	{
		GeneralMatrix* gm;
		int search(const BaseMatrix*) const;
	public:
		~ReturnMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		friend class BaseMatrix;
		ReturnMatrix(const ReturnMatrix& tm) : BaseMatrix(), gm(tm.gm) {}
		ReturnMatrix(const GeneralMatrix* gmx) : gm((GeneralMatrix*&)gmx) {}
		//   ReturnMatrix(GeneralMatrix&);
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(ReturnMatrix)
	};


	// ************************** submatrices ******************************/

	/// A submatrix of a matrix.
	/// \internal
	class LibGnut_LIBRARY_EXPORT GetSubMatrix : public NegatedMatrix
	{
		int row_skip;
		int row_number;
		int col_skip;
		int col_number;
		bool IsSym;

		GetSubMatrix
		(const BaseMatrix* bmx, int rs, int rn, int cs, int cn, bool is)
			: NegatedMatrix(bmx),
			row_skip(rs), row_number(rn), col_skip(cs), col_number(cn), IsSym(is) {}
		void SetUpLHS();
		friend class BaseMatrix;
	public:
		GetSubMatrix(const GetSubMatrix& g)
			: NegatedMatrix(g.bm), row_skip(g.row_skip), row_number(g.row_number),
			col_skip(g.col_skip), col_number(g.col_number), IsSym(g.IsSym) {}
		~GetSubMatrix() {}
		GeneralMatrix* Evaluate(MatrixType mt = MatrixTypeUnSp);
		void operator=(const BaseMatrix&);
		void operator+=(const BaseMatrix&);
		void operator-=(const BaseMatrix&);
		void operator=(const GetSubMatrix& m) { operator=((const BaseMatrix&)m); }
		void operator<<(const BaseMatrix&);
		void operator<<(const double*);                // copy from array
		void operator<<(const float*);                // copy from array
		void operator<<(const int*);                 // copy from array
		MatrixInput operator<<(double);                // for loading a list
		MatrixInput operator<<(float);                // for loading a list
		MatrixInput operator<<(int f);
		void operator=(Real);                        // copy from constant
		void operator+=(Real);                       // add constant
		void operator-=(Real r) { operator+=(-r); }  // subtract constant
		void operator*=(Real);                       // multiply by constant
		void operator/=(Real r) { operator*=(1.0 / r); } // divide by constant
		void inject(const GeneralMatrix&);           // copy stored els only
		void Inject(const GeneralMatrix& GM_A) { inject(GM_A); }
		MatrixBandWidth bandwidth() const;
		NEW_DELETE(GetSubMatrix)
	};

	// ******************** linear equation solving ****************************/

	/// A class for finding A.i() * B.
	/// This is supposed to choose the appropriate method depending on the
	/// type A. Not very satisfactory as it doesn't know about Cholesky for
	/// for positive definite matrices.
	class LibGnut_LIBRARY_EXPORT LinearEquationSolver : public BaseMatrix
	{
		GeneralMatrix* gm;
		int search(const BaseMatrix*) const { return 0; }
		friend class BaseMatrix;
	public:
		LinearEquationSolver(const BaseMatrix& bm);
		~LinearEquationSolver() { delete gm; }
		void cleanup() { delete gm; }
		GeneralMatrix* Evaluate(MatrixType) { return gm; }
		// probably should have an error message if MatrixType != UnSp
		NEW_DELETE(LinearEquationSolver)
	};

	// ************************** matrix input *******************************/

	/// Class for reading values into a (small) matrix within a program.
	/// \internal
	/// Is able to detect a mismatch in the number of elements.

	class LibGnut_LIBRARY_EXPORT MatrixInput
	{
		int n;                  // number values still to be read
		Real* r;                // pointer to next location to be read to
	public:
		MatrixInput(const MatrixInput& mi) : n(mi.n), r(mi.r) {}
		MatrixInput(int nx, Real* rx) : n(nx), r(rx) {}
		~MatrixInput();
		MatrixInput operator<<(double);
		MatrixInput operator<<(float);
		MatrixInput operator<<(int f);
		friend class GeneralMatrix;
	};



	// **************** a very simple integer array class ********************/

	/// A very simple integer array class.
	/// A minimal array class to imitate a C style array but giving dynamic storage
	/// mostly intended for internal use by newmat.
	/// Probably to be replaced by a templated class when I start using templates.

	class LibGnut_LIBRARY_EXPORT SimpleIntArray : public Janitor
	{
	protected:
		int* a;                    ///< pointer to the array
		int n;                     ///< length of the array
	public:
		SimpleIntArray(int xn);    ///< build an array length xn
		SimpleIntArray() : a(0), n(0) {}  ///< build an array length 0
		~SimpleIntArray();         ///< return the space to memory
		int& operator[](int i);    ///< access element of the array - start at 0
		int operator[](int i) const;
		///< access element of constant array
		void operator=(int ai);    ///< set the array equal to a constant
		void operator=(const SimpleIntArray& b);
		///< copy the elements of an array
		SimpleIntArray(const SimpleIntArray& b);
		///< make a new array equal to an existing one
		int Size() const { return n; }
		///< return the size of the array
		int size() const { return n; }
		///< return the size of the array
		int* Data() { return a; }  ///< pointer to the data
		const int* Data() const { return a; }  ///< pointer to the data
		int* data() { return a; }  ///< pointer to the data
		const int* data() const { return a; }  ///< pointer to the data
		const int* const_data() const { return a; }  ///< pointer to the data
		void resize(int i, bool keep = false);
		///< change length, keep data if keep = true
		void ReSize(int i, bool keep = false) { resize(i, keep); }
		///< change length, keep data if keep = true
		void resize_keep(int i) { resize(i, true); }
		///< change length, keep data
		void cleanup() { resize(0); }   ///< set length to zero
		void CleanUp() { resize(0); }   ///< set length to zero
		NEW_DELETE(SimpleIntArray)
	};

	// ********************** C subscript classes ****************************

	/// Let matrix simulate a C type two dimensional array
	class LibGnut_LIBRARY_EXPORT RealStarStar
	{
		Real** a;
	public:
		RealStarStar(Matrix& A);
		~RealStarStar() { delete[] a; }
		operator Real** () { return a; }
	};

	/// Let matrix simulate a C type const two dimensional array
	class LibGnut_LIBRARY_EXPORT ConstRealStarStar
	{
		const Real** a;
	public:
		ConstRealStarStar(const Matrix& A);
		~ConstRealStarStar() { delete[] a; }
		operator const Real** () { return a; }
	};

	// *************************** exceptions ********************************/

	/// Not positive definite exception.
	class LibGnut_LIBRARY_EXPORT NPDException : public Runtime_error
	{
	public:
		static unsigned long Select;
		NPDException(const GeneralMatrix&);
	};

	/// Covergence failure exception.
	class LibGnut_LIBRARY_EXPORT ConvergenceException : public Runtime_error
	{
	public:
		static unsigned long Select;
		ConvergenceException(const GeneralMatrix& A);
		ConvergenceException(const char* c);
	};

	/// Singular matrix exception.
	class LibGnut_LIBRARY_EXPORT SingularException : public Runtime_error
	{
	public:
		static unsigned long Select;
		SingularException(const GeneralMatrix& A);
	};

	/// Real overflow exception.
	class LibGnut_LIBRARY_EXPORT OverflowException : public Runtime_error
	{
	public:
		static unsigned long Select;
		OverflowException(const char* c);
	};

	/// Miscellaneous exception (details in character string). 
	class LibGnut_LIBRARY_EXPORT ProgramException : public Logic_error
	{
	protected:
		ProgramException();
	public:
		static unsigned long Select;
		ProgramException(const char* c);
		ProgramException(const char* c, const GeneralMatrix&);
		ProgramException(const char* c, const GeneralMatrix&, const GeneralMatrix&);
		ProgramException(const char* c, MatrixType, MatrixType);
	};

	/// Index exception.
	class LibGnut_LIBRARY_EXPORT IndexException : public Logic_error
	{
	public:
		static unsigned long Select;
		IndexException(int i, const GeneralMatrix& A);
		IndexException(int i, int j, const GeneralMatrix& A);
		// next two are for access via element function
		IndexException(int i, const GeneralMatrix& A, bool);
		IndexException(int i, int j, const GeneralMatrix& A, bool);
	};

	/// Cannot convert to vector exception.
	class LibGnut_LIBRARY_EXPORT VectorException : public Logic_error
	{
	public:
		static unsigned long Select;
		VectorException();
		VectorException(const GeneralMatrix& A);
	};

	/// A matrix is not square exception.
	class LibGnut_LIBRARY_EXPORT NotSquareException : public Logic_error
	{
	public:
		static unsigned long Select;
		NotSquareException(const GeneralMatrix& A);
		NotSquareException();
	};

	/// Submatrix dimension exception.
	class LibGnut_LIBRARY_EXPORT SubMatrixDimensionException : public Logic_error
	{
	public:
		static unsigned long Select;
		SubMatrixDimensionException();
	};

	/// Incompatible dimensions exception.
	class LibGnut_LIBRARY_EXPORT IncompatibleDimensionsException : public Logic_error
	{
	public:
		static unsigned long Select;
		IncompatibleDimensionsException();
		IncompatibleDimensionsException(const GeneralMatrix&);
		IncompatibleDimensionsException(const GeneralMatrix&, const GeneralMatrix&);
	};

	/// Not defined exception.
	class LibGnut_LIBRARY_EXPORT NotDefinedException : public Logic_error
	{
	public:
		static unsigned long Select;
		NotDefinedException(const char* op, const char* matrix);
	};

	/// Cannot build matrix with these properties exception.
	class LibGnut_LIBRARY_EXPORT CannotBuildException : public Logic_error
	{
	public:
		static unsigned long Select;
		CannotBuildException(const char* matrix);
	};


	/// Internal newmat exception - shouldn't happen.
	class LibGnut_LIBRARY_EXPORT InternalException : public Logic_error
	{
	public:
		static unsigned long Select;          // for identifying exception
		InternalException(const char* c);
	};

	// ************************ functions ************************************ //

	LibGnut_LIBRARY_EXPORT bool operator==(const GeneralMatrix& A, const GeneralMatrix& B);
	LibGnut_LIBRARY_EXPORT bool operator==(const BaseMatrix& A, const BaseMatrix& B);
	LibGnut_LIBRARY_EXPORT inline bool operator!=(const GeneralMatrix& A, const GeneralMatrix& B)
	{
		return !(A == B);
	}
	LibGnut_LIBRARY_EXPORT inline bool operator!=(const BaseMatrix& A, const BaseMatrix& B)
	{
		return !(A == B);
	}

	// inequality operators are dummies included for compatibility
	// with STL. They throw an exception if actually called.
	LibGnut_LIBRARY_EXPORT inline bool operator<=(const BaseMatrix& A, const BaseMatrix&)
	{
		A.IEQND(); return true;
	}
	LibGnut_LIBRARY_EXPORT inline bool operator>=(const BaseMatrix& A, const BaseMatrix&)
	{
		A.IEQND(); return true;
	}
	LibGnut_LIBRARY_EXPORT inline bool operator<(const BaseMatrix& A, const BaseMatrix&)
	{
		A.IEQND(); return true;
	}
	LibGnut_LIBRARY_EXPORT inline bool operator>(const BaseMatrix& A, const BaseMatrix&)
	{
		A.IEQND(); return true;
	}

	LibGnut_LIBRARY_EXPORT bool is_zero(const BaseMatrix& A);
	LibGnut_LIBRARY_EXPORT inline bool IsZero(const BaseMatrix& A) { return is_zero(A); }

	LibGnut_LIBRARY_EXPORT Real dotproduct(const Matrix& A, const Matrix& B);
	LibGnut_LIBRARY_EXPORT Matrix crossproduct(const Matrix& A, const Matrix& B);
	LibGnut_LIBRARY_EXPORT ReturnMatrix crossproduct_rows(const Matrix& A, const Matrix& B);
	LibGnut_LIBRARY_EXPORT ReturnMatrix crossproduct_columns(const Matrix& A, const Matrix& B);

	LibGnut_LIBRARY_EXPORT inline Real DotProduct(const Matrix& A, const Matrix& B)
	{
		return dotproduct(A, B);
	}
	LibGnut_LIBRARY_EXPORT inline Matrix CrossProduct(const Matrix& A, const Matrix& B)
	{
		return crossproduct(A, B);
	}
	LibGnut_LIBRARY_EXPORT inline ReturnMatrix CrossProductRows(const Matrix& A, const Matrix& B)
	{
		return crossproduct_rows(A, B);
	}
	LibGnut_LIBRARY_EXPORT inline ReturnMatrix CrossProductColumns(const Matrix& A, const Matrix& B)
	{
		return crossproduct_columns(A, B);
	}

	LibGnut_LIBRARY_EXPORT void newmat_block_copy(int n, Real* from, Real* to);

	// ********************* friend functions ******************************** //

	// Functions declared as friends - G++ wants them declared externally as well

	LibGnut_LIBRARY_EXPORT bool Rectangular(MatrixType a, MatrixType b, MatrixType c);
	LibGnut_LIBRARY_EXPORT bool Compare(const MatrixType&, MatrixType&);
	LibGnut_LIBRARY_EXPORT Real dotproduct(const Matrix& A, const Matrix& B);
	LibGnut_LIBRARY_EXPORT SPMatrix SP(const BaseMatrix&, const BaseMatrix&);
	LibGnut_LIBRARY_EXPORT KPMatrix KP(const BaseMatrix&, const BaseMatrix&);
	LibGnut_LIBRARY_EXPORT ShiftedMatrix operator+(Real f, const BaseMatrix& BM);
	LibGnut_LIBRARY_EXPORT NegShiftedMatrix operator-(Real, const BaseMatrix&);
	LibGnut_LIBRARY_EXPORT ScaledMatrix operator*(Real f, const BaseMatrix& BM);

	// ********************* inline functions ******************************** //

	LibGnut_LIBRARY_EXPORT inline LogAndSign log_determinant(const BaseMatrix& B)
	{
		return B.log_determinant();
	}
	LibGnut_LIBRARY_EXPORT inline LogAndSign LogDeterminant(const BaseMatrix& B)
	{
		return B.log_determinant();
	}
	LibGnut_LIBRARY_EXPORT inline Real determinant(const BaseMatrix& B)
	{
		return B.determinant();
	}
	LibGnut_LIBRARY_EXPORT inline Real Determinant(const BaseMatrix& B)
	{
		return B.determinant();
	}
	LibGnut_LIBRARY_EXPORT inline Real sum_square(const BaseMatrix& B) { return B.sum_square(); }
	LibGnut_LIBRARY_EXPORT inline Real SumSquare(const BaseMatrix& B) { return B.sum_square(); }
	LibGnut_LIBRARY_EXPORT inline Real norm_Frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
	LibGnut_LIBRARY_EXPORT inline Real norm_frobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
	LibGnut_LIBRARY_EXPORT inline Real NormFrobenius(const BaseMatrix& B) { return B.norm_Frobenius(); }
	LibGnut_LIBRARY_EXPORT inline Real trace(const BaseMatrix& B) { return B.trace(); }
	LibGnut_LIBRARY_EXPORT inline Real Trace(const BaseMatrix& B) { return B.trace(); }
	LibGnut_LIBRARY_EXPORT inline Real sum_absolute_value(const BaseMatrix& B)
	{
		return B.sum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline Real SumAbsoluteValue(const BaseMatrix& B)
	{
		return B.sum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline Real sum(const BaseMatrix& B)
	{
		return B.sum();
	}
	LibGnut_LIBRARY_EXPORT inline Real Sum(const BaseMatrix& B)
	{
		return B.sum();
	}
	LibGnut_LIBRARY_EXPORT inline Real maximum_absolute_value(const BaseMatrix& B)
	{
		return B.maximum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline Real MaximumAbsoluteValue(const BaseMatrix& B)
	{
		return B.maximum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline Real minimum_absolute_value(const BaseMatrix& B)
	{
		return B.minimum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline Real MinimumAbsoluteValue(const BaseMatrix& B)
	{
		return B.minimum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline Real maximum(const BaseMatrix& B) { return B.maximum(); }
	LibGnut_LIBRARY_EXPORT inline Real Maximum(const BaseMatrix& B) { return B.maximum(); }
	LibGnut_LIBRARY_EXPORT inline Real minimum(const BaseMatrix& B) { return B.minimum(); }
	LibGnut_LIBRARY_EXPORT inline Real Minimum(const BaseMatrix& B) { return B.minimum(); }
	LibGnut_LIBRARY_EXPORT inline Real norm1(const BaseMatrix& B) { return B.norm1(); }
	LibGnut_LIBRARY_EXPORT inline Real Norm1(const BaseMatrix& B) { return B.norm1(); }
	LibGnut_LIBRARY_EXPORT inline Real norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
	LibGnut_LIBRARY_EXPORT inline Real Norm1(RowVector& RV) { return RV.maximum_absolute_value(); }
	LibGnut_LIBRARY_EXPORT inline Real norm_infinity(const BaseMatrix& B) { return B.norm_infinity(); }
	LibGnut_LIBRARY_EXPORT inline Real NormInfinity(const BaseMatrix& B) { return B.norm_infinity(); }
	LibGnut_LIBRARY_EXPORT inline Real norm_infinity(ColumnVector& CV)
	{
		return CV.maximum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline Real NormInfinity(ColumnVector& CV)
	{
		return CV.maximum_absolute_value();
	}
	LibGnut_LIBRARY_EXPORT inline bool IsZero(const GeneralMatrix& A) { return A.IsZero(); }
	LibGnut_LIBRARY_EXPORT inline bool is_zero(const GeneralMatrix& A) { return A.is_zero(); }


	LibGnut_LIBRARY_EXPORT inline MatrixInput MatrixInput::operator<<(int f) { return *this << (Real)f; }
	LibGnut_LIBRARY_EXPORT inline MatrixInput GeneralMatrix::operator<<(int f) { return *this << (Real)f; }
	LibGnut_LIBRARY_EXPORT inline MatrixInput BandMatrix::operator<<(int f) { return *this << (Real)f; }
	LibGnut_LIBRARY_EXPORT inline MatrixInput GetSubMatrix::operator<<(int f) { return *this << (Real)f; }
	LibGnut_LIBRARY_EXPORT inline ReversedMatrix BaseMatrix::Reverse() const { return reverse(); }
	LibGnut_LIBRARY_EXPORT inline RowedMatrix BaseMatrix::AsRow() const { return as_row(); }
	LibGnut_LIBRARY_EXPORT inline ColedMatrix BaseMatrix::AsColumn() const { return as_column(); }
	LibGnut_LIBRARY_EXPORT inline DiagedMatrix BaseMatrix::AsDiagonal() const { return as_diagonal(); }
	LibGnut_LIBRARY_EXPORT inline MatedMatrix BaseMatrix::AsMatrix(int m, int n) const
	{
		return as_matrix(m, n);
	}
	LibGnut_LIBRARY_EXPORT inline GetSubMatrix BaseMatrix::SubMatrix(int fr, int lr, int fc, int lc) const
	{
		return submatrix(fr, lr, fc, lc);
	}
	LibGnut_LIBRARY_EXPORT inline GetSubMatrix BaseMatrix::SymSubMatrix(int f, int l) const
	{
		return sym_submatrix(f, l);
	}
	LibGnut_LIBRARY_EXPORT inline GetSubMatrix BaseMatrix::Row(int f) const { return row(f); }
	LibGnut_LIBRARY_EXPORT inline GetSubMatrix BaseMatrix::Rows(int f, int l) const { return rows(f, l); }
	LibGnut_LIBRARY_EXPORT inline GetSubMatrix BaseMatrix::Column(int f) const { return column(f); }
	LibGnut_LIBRARY_EXPORT inline GetSubMatrix BaseMatrix::Columns(int f, int l) const
	{
		return columns(f, l);
	}
	LibGnut_LIBRARY_EXPORT inline Real BaseMatrix::AsScalar() const { return as_scalar(); }

	LibGnut_LIBRARY_EXPORT inline ReturnMatrix GeneralMatrix::ForReturn() const { return for_return(); }

	LibGnut_LIBRARY_EXPORT inline void swap(Matrix& A, Matrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(SquareMatrix& A, SquareMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(nricMatrix& A, nricMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(UpperTriangularMatrix& A, UpperTriangularMatrix& B)
	{
		A.swap(B);
	}
	LibGnut_LIBRARY_EXPORT inline void swap(LowerTriangularMatrix& A, LowerTriangularMatrix& B)
	{
		A.swap(B);
	}
	LibGnut_LIBRARY_EXPORT inline void swap(SymmetricMatrix& A, SymmetricMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(DiagonalMatrix& A, DiagonalMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(RowVector& A, RowVector& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(ColumnVector& A, ColumnVector& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(CroutMatrix& A, CroutMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(BandMatrix& A, BandMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(UpperBandMatrix& A, UpperBandMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(LowerBandMatrix& A, LowerBandMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(SymmetricBandMatrix& A, SymmetricBandMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(BandLUMatrix& A, BandLUMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(IdentityMatrix& A, IdentityMatrix& B) { A.swap(B); }
	LibGnut_LIBRARY_EXPORT inline void swap(GenericMatrix& A, GenericMatrix& B) { A.swap(B); }

#ifdef OPT_COMPATIBLE                    // for compatibility with opt++

	inline Real Norm2(const ColumnVector& CV) { return CV.norm_Frobenius(); }
	inline Real Dot(ColumnVector& CV1, ColumnVector& CV2)
	{
		return dotproduct(CV1, CV2);
	}

#endif


#ifdef use_namespace
}
#endif


#endif

// body file: newmat1.cpp
// body file: newmat2.cpp
// body file: newmat3.cpp
// body file: newmat4.cpp
// body file: newmat5.cpp
// body file: newmat6.cpp
// body file: newmat7.cpp
// body file: newmat8.cpp
// body file: newmatex.cpp
// body file: bandmat.cpp
// body file: submat.cpp



///@}




